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Abstract. We present here both analytical and numerical results of hydrodynamic stability investigations of rota- 
tionally supported circumstellar flows using the shearing box formalism. Asymptotic scaling arguments justifying 
the shearing box approximation are systematically derived, showing that there exist two limits which we call 
small shearing box (SSB) and large shearing box (LSB). The physical meaning of these two limits and their rela- 
tionship to model equations implemented by previous investigators are discussed briefly. Two dimensional (2D) 
dynamics of the SSB are explored and shown to contain transiently growing (TG) linear modes, whose nature is 
first discussed within the context of linear theory. The fully nonlinear regime in 2D is investigated numerically 
for very high Reynolds (Re) numbers. Solutions exhibiting long-term dynamical activity are found and manifest 
episodic but recurrent TG behavior and these are associated with the formation and long-term survival of coherent 
vortices. The life-time of this spatio-temporal complexity depends on the Re number and the strength and nature 
of the initial disturbance. The dynamical activity in finite Re solutions ultimately decays with a characteristic 
time increasing with Re. However, for large enough Re and appropriate initial perturbation, a large number of 
TG episodes recur before any viscous decay begins to clearly manifest itself. In cases where Re = oo nominally 
(i.e. any dissipation resulting only from numerical truncation errors), the dynamical activity persists for the entire 
duration of the simulation (hundreds of box orbits). Because the SSB approximation used here is equivalent to 
a 2D incompressible flow, the dynamics can not depend on the Coriolis force. Therefore, three dimensional (3D) 
simulations are needed in order to decide if this force indeed suppresses nonlinear hydrodynamical instability in 
rotationally supported disks in the shearing box approximation, and if recurrent TG behavior can still persist in 
three dimensions as well - possibly giving rise to a subcritical transition to long-term spatio-temporal complexity. 
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1. Introduction 

Accretion of matter, endowed with a significant amount 
of angular momentum, onto relatively compact objects, 
is thought to occur in a variety of important astrophys- 
ical systems, including close binary stars, young stellar 
objects and active galactic nuclei. If the accreting matter 
is able to cool efficiently enough, the resulting configura- 
tion of the flow is a vertically thin, almost-Keplerian (that 
is, essentially rotationally supported), accretion disk. The 
accretion process (transport of mass inward) can occur 
only if some dissipative mechanism is present, giving rise 
to transport of angular momentum outward. Since the mi- 
croscopic viscosity of the fluids in question is very small, 
and thus the Re numbers of such flows are literally as- 
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tronomical, anomalous transport of some kind must be 
invoked to account for the time scales and magnitude of 
the energy output of the accreting sources. The pioneer- 
ing works of Shakura & Sunyaev (1973) and Lyndcn-Bcll 
& Pringle (1972) suggested that the needed anomalous 
transport can be naturally achieved if the accretion flow 
is turbulent and, moreover, proposed to circumvent the 
complexities and our lack of understanding of turbulence 
by introducing the now famous a model of accretion disks 
(in the language of turbulence modeling it is "a zero- 
equation, one-parameter model"). This simple idea has 
been extremely fruitful, giving rise to a large number of 
successful interpretations of observational results (see e.g. 
Pringle 1981 and Frank, King & Raine 2002 for reviews). 

The theoretical question of the physical origin of 
accretion disk turbulence has however remained essen- 
tially unanswered until Balbus & Hawley (1991) dis- 
covered that weak magnetic fields destabilize differential 
Keplerian rotation, i.e. the magneto-rotational instabil- 
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ity (Velikhov 1959, Chandrasekhar 1960), hereafter MRI, 
can operate in rotationally supported magnetized accre- 
tion disks, rendering them linearly unstable. This insta- 
bility was subsequently recognized as the source of accre- 
tion disk magneto-hydrodynamic (MHD) turbulence, giv- 
ing rise to the needed angular momentum transport and 
also to dynamo action (Hawley, Gammie & Balbus 1995; 
Brandenburg et al. 1995). For a recent review of this topic 
see Balbus (2003). The discovery by Balbus & Hawley is 
very important and significant because it demonstrates 
the existence of linear instability in accretion disk mod- 
els with angular velocity profiles that are otherwise stable 
according to the Rayleigh criterion. Although some ideas 
have been put forward in this context, e.g. the possibility 
of an instability resulting from the angular velocity be- 
ing nonconstant on cylinders (Knobloch & Spruit 1986 
Kluzniak & Kita 1997, Regev & Gitelman 2002, Urpin 
2003) or more general baroclinic instabilities (Klahr & 
Bodcnheimer 2003), no hydrodynamic instability of any 
kind has ever been explicitly shown to exist. In the absence 
of a known linear hydrodynamic instability, the turbulent 
viscosity needed for angular momentum transport in ac- 
cretion disks has acquired in the community a somewhat 
"mysterious" character, similar to other ignorance-driven 
concepts in the history of astrophysics (and physics). Thus 
the MRI has been generally accepted with considerable 
enthusiasm and, as it often happens in cases where a suc- 
cessful idea is advanced, raised to the level of an exclu- 
sive paradigm. Moreover, the very possibility of the occur- 
rence of hydrodynamic turbulence (resulting from linear 
or even nonlinear hydrodynamical instabilities) has been 
questioned by Balbus, Hawley & Stone (1996) (see also the 
detailed study of Hawley, Balbus & Winters 1999), here- 
after BHSW, on the basis of the results of fully nonlinear 
finite-difference numerical simulations in the 3D shearing 
box approximation (see below). 

The extreme displeasure felt by most astrophysicists 
with the absence of linear hydrodynamical disk instabil- 
ities, and their reluctance in accepting that disks may 
still be hydrodynamically turbulent, seems to be detached 
from what has been known for the past century on hydro- 
dynamical stability of shear flows. It is well known (see e.g. 
Joseph 1972, Drazin & Reid 1983) that the linear stabil- 
ity properties of essentially all the canonical shear flows, 
rotating or not, do not reproduce well their laboratory 
behavior. In almost all such flows transition to turbulence 
is subcritical, that is, it occurs at Reynolds numbers (Re) 
that are significantly lower (and also dependent on the 
perturbation strength) than the critical values for linear 
instability Re c . In some instances, like the pipe Poiseulle, 
plane Couette or Taylor-Couette (in the narrow gap limit) 
flows, this problem is particularly acute as linear analysis 
of perturbations on these basic flow profiles predicts sta- 
bility (even decay) of infinitesimal perturbations, that is, 
formally Re c = oo. Laboratory experiments demonstrate, 
however, full-on transition to turbulence via some insta- 
bility mechanism that is neither fully understood nor even 



uniquely identified. Thus various ideas in this context are 
generally referred to as transition scenarios (see below) . 

On the theoretical side, it has been known since the 
work of Orr (1907) that even though a basic shear flow 
may be linearly stable, some perturbations (i.e. initial 
conditions) may exhibit significant transient growth (TG) 
within the linear regime, before ultimately decaying. This 
fact prompted a number of researchers to examine the 
possibility of subcritical transition to turbulence, with the 
linearly stable flow finding a way to bypass the usual (i.e. 
via linear instability) route to turbulence. We are unable 
to mention here all the contributions on this subject (in 
particular, those published before the early 1990s) and re- 
fer the interested reader to one of the reviews cited below, 
but only note that an application of TG to an astrophys- 
ical problem has already been implemented by Goldreich 
& Lynden-Bell (1965). 

A significant step forward in the advancement of the 
bypass transistion idea has been made in the early 1990s. 
It came with the clarification that TG, and thus the pos- 
sibility of a bypass transition, is the result of the lin- 
earized system's operator properties (Butler & Farrell 
1992, Reddy & Henningson 1993). The essence of this 
idea lies in the observation that linear stability analysis of 
shear flows gives rise to a nonnormal operator(i.e. an op- 
erator which does not commute with its adjoint), that is, 
one whose eigenfunctions are not necessarily orthogonal. 
This rather elementary mathematical fact may have pro- 
found consequences for hydrodynamic stability, because 
it explains how TG may naturally occur in linearly stable 
basic shear flows (see e.g. Trefethen et al. 1993). This ap- 
proach is frequently referred to as nonmodal because the 
linear stability analysis required to find TG cannot pro- 
ceed via the usual modal analysis, involving a boundary 
value problem and must consist of the examination of the 
original initial value problem. If this TG is strong enough 
then nonlinear interactions may come into play. This then 
could further lead to some kind of " nonlinear instability" 
creating a sustained complex dynamical state in which 
TG events persistently recur and resulting in a situation 
resembling turbulence. We shall refer to this phenomenon 
as recurrent transient growth (RTG) . The recent review by 
Grossmann (2000) contains an excellent account of these 
ideas. We refer the reader also to the comprehensive mod- 
ern book on shear flow instabilities (Schmid & Henningson 
2001) and relevant references within, where a detailed ac- 
count (a full chapter) on the different ideas on transition 
to turbulence in shear flows are detailed. 

The above developments recently have been brought to 
the attention of the astrophysical community by Ioannaou 
& Kakouris (2001) and Chagelishvili et al. (2003) in the 
context of accretion disks. These two studies are quite 
different in their approach - the former examines global 
behavior of disturbances in an accretion disk and the lat- 
ter is a local analysis in the framework of the shearing box 
(SB), also sometimes called the shearing sheet, approxima- 
tion - but they both utilize the basic idea mentioned above. 
It appears that although these works have not explicitly 
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demonstrated that a turbulent state of an accretion disk 
actually arises by purely hydrodynamical processes, they 
have convincingly suggested that such a process is feasible. 
The purpose of this (and subsequent) work is to carry this 
idea further and examine in detail the nonlinear regime in 
various limits. 

There is little doubt that whenever MRI occurs in a 
disk it leads to a turbulent flow with enhanced angular 
momentum transport outward. However, BHSW (see also 
Balbus 2003) put forth the claim that no purely hydro- 
dynamical process can destabilize a Keplerian accretion 
disk, at least within the SB approximation. This is based 
on fully nonlinear simulations (performed by BHSW) in 
which no (nonlinear) hydrodynamic instability was shown 
to manifest itself as long as the angular velocity profile 
was linearly stable according to the Rayleigh criterion. 
Furthermore, the same simulations of BHSW also show 
that the total kinetic energy of plane Couette flow and ro- 
tating flows which are unstable according to the Rayleigh 
criterion, grow essentially monotonically (after a short 
transient). By contrast, in Rayleigh-stable rotation pro- 
files the kinetic energy was observed to decay over the 
duration of the runs (typically 5-10 box orbit times). 

Since in the SB approximation it is the Coriolis term 
that actually distinguishes between simple plane Couette 
flow and the flow within a small segment of a rotating 
disk (see below), it is only natural that this force has been 
identified by BHSW as the stabilizing agent. 

These numerical results certainly constitute a serious 
challenge to the possibility of hydrodynamic turbulence in 
rotationally supported disks. On the other hand, we can 
state at least five observations which suggest that this 
matter is far from being settled. 

1. In some astrophysical systems containing such disks 
the conditions are such that MRI driven MHD turbu- 
lence is probably impossible. For example, accretion 
disks in dwarf novae during quiescence (Menou 2000) 
and around forming stars (Blaes & Balbus 1994, Sano 
et al. 2000) seem to be too resistive to support MHD 
turbulence. 

2. Recent laboratory Couette- Taylor experiments in the 
narrow gap limit, with linearly stable rotational angu- 
lar velocity profiles (like in Keplerian disks) indicate 
the development of turbulence in such flows (Richard 
2001, see also Richard & Zahn 1999). 

3. In a recent thoughtful examination of the problem, 
Longaretti (2002) proposes that the above conclusion 
reached by BHSW, about the role of the Coriolis force 
in suppressing the subcritical transition (observed in 
non-rotating plane Couette flows in the lab and in 
BHSW simulations), may be premature. He speculates 
that their findings may stem from the inability (of the 
simulations) to resolve some of the essential dynamics 
within these rotating flows (see however Balbus 2003). 

4. There exist numerical simulations of plane Couette 
flows, which have revealed the details of what can 
be called "subcritical turbulence" in these flows 



(Schmiegel & Eckhardt, 1997; Eckhardt et al. 1998). 
These calculations, which are the only ones (as far 
as we know) of this kind and extend for reasonably 
long integration times, indicate that above a critical 
Reynolds number there appears dynamical activity, re- 
sulting probably from repeated TG events, brought 
about by nonlinear interaction and persisting much 
longer than the viscous decay time, i.e. RTG. This 
nature of the subcritical bypass transition appears to 
be very different from a "usual" supercritical transi- 
tion to turbulence via linear instability. In particular 
the outcome critically depends on the initial pertur- 
bation, in addition to the Reynolds number. Similar 
dependence is observed in experiments (see Grossman 
2000). It is therefore quite surprising that in BHSW 
the non-rotating Couette flow (which is linearly stable 
and hence must have a subcritical transition and thus 
probably RTG) and the rotating Rayleigh (linearly) 
unstable flows exhibit a similar perturbation energy 
growth, at least during the relatively short duration of 
the BHSW simulations and for a particular choice of 
the initial perturbation. 
5. In relation to the last point, the simulations of BHSW 
only consider one type of initial conditions to seed 
their flows and they follow its evolution for, at most, 
a dozen orbit times. Perturbation spectra are surely 
myriad and it is possible that conditions leading to a 
subcritical transition into turbulence may have been 
overlooked by BHSW. 

Thus, it seems that continued study of purely hydrody- 
namical processes in disks still remains viable and worth- 
while (in the words of Balbus, Hawley & Stone 1996) and, 
in particular, the issue of hydrodynamical stability war- 
rants perhaps a new look. 

In this paper we will report on our first effort in this di- 
rection. Our goal will be to critically evaluate, using high 
resolution and long duration (hundreds of box orbit times, 
that is, Kepler rotational periods at the radius of the box 
location) numerical experiments in the nonlinear regime, 
whether the linear TG mechanism is sufficient to induce 
a bypass transition into a persistently dynamically active 
state, like the one described in item (4) above. We shall 
employ spectral methods (wherever possible) for our sim- 
ulations, as these methods are generally regarded in the 
fluid-dynamical community as being more accurate and 
reliable than finite difference schemes. For the time being 
we shall focus just on the nature of the instability and the 
possibility of a sustained spatio-temporal complexity. The 
issue whether in such a state there is an enhanced angular 
momentum transport outward (as it must be in accretion 
disks) is an important one (see Cabot 1966, Balbus 2003), 
but its exploration must await a conclusive result on the 
instability itself and thus we shall not deal with it in the 
present work. 

This paper is organized in the following way. In Section 
2 we examine and discuss the nondimensional SB equa- 
tions, as appropriate for the case studies here, in which 
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the size of the box is much smaller than the disk's vertical 
scale height. This is one of the limits of the SB approxima- 
tion, whose systematic and general derivation is detailed 
in Appendix A. Section 3 deals with linear theory and in 
particular we apply there the concept of sheared coordi- 
nates, with whose help the results of Chagelishvili et al. 
(2003) (the TG of linear modes) can be elegantly recov- 
ered. In section 4 we report the results of nonlinear 2D 
numerical simulation in this limit of the SB approxima- 
tion. Although it has been recognized that nonmodal TG 
may give rise to very large amplification in 3D, also 2D 
perturbations can grow quite significantly (Trefethen et 
al. 1993). In two-dimensional viscous incompressible flows 
any dynamical activity must ultimately decay (see section 
4.2 below), however for large enough Re numbers we ex- 
pect the decay times to be very long (since they should 
behave as oc Re 2/3 , Yecko 2004) so as to allow for the TG 
episodes to recur many times. 

Thus, in this initial study we limit ourselves to 2D 
and we postpone to subsequent investigations (which are 
already in progress) 3D simulations of this and the other 
limit of the SB approximation. 

2. The Shearing Box Approximation Equations: 

Numerical simulations of rotationally supported flows at 
very large Reynolds numbers must focus on only very 
small disk sections if high spatial resolution is imperative. 
This is clearly the case when one looks for possible insta- 
bilities and therefore most such studies have been done 
within the shearing box (or sheet) approximation. The 
essence of this approximation is local in approach, that 
is, the equations are valid in a small region (a Cartesian 
box) about a typical point in the disk. In such a box, 
a steady flow consisting of a linear shear velocity profile 
solves the equations and one can consider it as a basic 
flow and perturb around it. The equations for the per- 
turbations (linearized or not) are then referred to as the 
shearing box (SB) approximation equations. Most numeri- 
cal simulations of Balbus, Hawley and collaborators in the 
MRI context have been done within this formalism as have 
the recent studies, mentioned above, on the feasibility of a 
hydrodynamical instability. Although the SB approxima- 
tion is not new (Goldreich & Lynden Bell 1965) we have 
not found in the literature a systematic derivation of this 
approximation. Thus, in the purpose of resolving the con- 
fusion (of the present authors and possibly also others) we 
give here in Appendix A such a derivation, giving rise to 
two limits of this approximation. 

In the first one, a box whose size is much smaller 
than the smallest scale height in the disk (usually the one 
in the vertical direction) is considered, and thus the un- 
perturbed state of linear shear may be considered homo- 
geneous. The perturbations are then incompressible and 
acoustics are ruled out. Most previously relevant work 
have been done in this approximation, which we call small 
shearing box (SSB). For the sake of completeness we re- 
peat here the non-dimcnsionalized SSB equations derived 



in Appendix A. In writing the equations below, we have 
deviated slightly from the notation presented in Appendix 
A, 



V • u = 0, 

dtu + 2Axd v u — 2H,qv + u • Vu = — 



d t v + 2Axd y v + 2(tt + A)u + u • Vw = - 



2Axd y w 



d t w 

d t p' + 2Axd v p' 



+ u • Vw = 
+ u • Vp' = 0, 



d x p 
l + p" 

d yP 
1+p" 
d z p + zp' 
l + o' 



(1) 
(2) 

(3) 

(4) 
(5) 

and the identifications to the terms in Appendix A are: 
u = u', u = u' x ,v = u' u ,w = u' z . We also specifically 
assume that the disk is exactly Keplerian which means 
that fio = 1 and 2A — —3/2. It also means, by l|A.20|) . 
that d z Pb = —z and d x Pb — 0. Additionally, even though 
fio = 1; we retain this symbol throughout the calculation 
for the purpose of flagging the Coriolis terms. 

In this paper we exclusively use the SSB approxima- 
tion and assume that the initial density disturbance is 
everywhere zero, i.e. p'(t = 0) = 0. From © it follows 
that p' = for all subsequent times. Thus, © becomes 
irrelevant and from here on out it shall not be referenced. 

The large shearing box (whose size is of the order the 
disk thickness) equations are also given in Appendix A 
and they will be treated in our future works. 

3. SSB linear theory - a review 

For the sake of economy we reintroduce the parameter q 
(with £lo = 1) as defined by IjA.llfl to replace — 2A in the 
linearized expressions of 111 II) : 



0, 



(A 



d x u + d y v + d z w 

(dt — qxd y )u — 2Q.qv 

qxd y )v + (2Qq — q)u 

(d t - qxd y )w 

Since this flow is exactly Keplerian, q ■■ 
of reference: for solid body rotation q = 
constant specific angular momentum q 



-d x p, 
-d v p, 
-d z p. 



(0) 
(7) 
(8) 
(9) 

= 3/2. As a point 
0, for flows with 
= 2. These linear 



equations describe simple incompressible flow with a lin- 
ear shear profile in a rotating frame: that is, linearized 
rotating plane Couette flow (Nagata, 1986). 

For two-dimensional disturbances (i.e. when d z = 0, 
w = 0), the equations governing the linear evolution sim- 
plify. With the vertical component of vorticity defined by, 

£ = d x v - d y u, (10) 

it follows from (JTJ and (|SJ that it is conserved along the 
shear flow, 

(dt ~ qxd y )t = 0. (11) 

Because of the incompressibihty of the 2D flow, the veloc- 
ity components may be derived from a stream function, 
"0, defined by 

u = -dy-ip, v = d x ip, (12) 
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and related to the vorticity via 



(13) 



We note also that the equations governing 2D dynam- 
ics are insensitive to the background rotation state even 
though the kinematics are sensitive to it. If the bound- 
ary conditions are independent of rotation too (as is the 
case here), the velocity distribution (that is the dynamics) 
cannot depend on the Coriolis term. This is a well known 
property of the governing equations and it is entirely due 
to the incompressibility of the flow which facilitates the 
streamfunction-vorticity formulation (see e.g. Batchelor 
1967, p. 178). The resulting equation set, even the non- 
linear one, becomes mathematically equivalent to the set 
appropriate to plane Couette flow. The results of all our 
2D calculations are based on solving equations l|ll|) and 
<|13l) with periodic boundary conditions (in the sheared 
frame, see below) and thus they do not depend in any 
way on the Coriolis force. 

To avoid any confusion with language, in this work 
we interchangeably refer to points with constant values 
of x and varying values of y as being in the streamwise 
direction and we refer to points with constant values of y 
and varying values of a; as being in the shearwise direction. 

Note that for three-dimensional disturbances one can 
obtain a single decoupled equation for the radial velocity 
component u after manipulation of I|6I9|) . 

(d t - qxdyf{dl + d 2 + d 2 z )u + 2flo(2r!o - q)d 2 u = 0. (14) 

Finally we point out that the same manipulations show 
that the pressure is related to the radial velocity by 



dip + d 2 z p = (d t - qxd y )d x u - (2f2 - q)d y u. 



(15) 



In contrast to 2D dynamics, 3D dynamics feel the effect of 
rotation through the explicit presence of the Coriolis term 
2fi (2^o - q) in lfl4) . 

We continue the analysis by reformulating the linear 
problem in terms of shearing coordinates ( SC for short) 
in the same way as implemented by Goldreich & Lynden- 
Bell (1965). To make sure our terminology is clear, we shall 
refer to the usual (untransformed) coordinate formulation 
of this problem to be in the non-shearing coordinates (NSC 
for short). The NSC is what an observer would see and, as 
such, we interchangeably refer to it as the observer frame. 

3.1. Shearing coordinates 

The coordinates are written into a frame which is shearing 
exactly as the background flow itself. This new coordinate 
system is formally defined as, 

X = x, Y = y + q(t - t )x, Z = z, T = t-t , (16) 

where t is some arbitrary reference time. Derivative op- 
erators are replaced in the following sense, 

d t = d T +qXd Y1 d x = d x + qTd Y , d y = d Y , d z — <9 Z .(17) 



The 3D equations that result from this transformation 
applied to equations (|6I9|I are, 



(d x + qTd Y ) u + d Y v + d z w = 0, 

d T u - 2f2 « = - (d x - 

d T v + (2Q - q)u = -d Y p, 

d T w = —d z p. 



qTd Y )p, 



(18a) 
(18b) 
(18c) 
(18d) 



The advantage in going over into this reference frame 
is that we explicitly remove the shear expression (which is 
proportional to x) from the governing equations. In return 
we receive a term that is proportional to time, T, in the 
resulting initial value problem (IVP). Rewriting (|14J) and 
JT2J in the SC gives, 



and 



(d x +qTd Y y 



d 2 Y+d 2 z 



q)dlu = 0, 



(19) 



(d 2 Y +d 2 z )p=[d T (d x +qTd Y )u-(2n Q -q)}u. (20) 

We remind the reader that this coordinate transforma- 
tion is volume preserving because its Jacobian is exactly 
one. 

We investigate the linear behavior in terms of the 
Fourier components of the disturbances. This is a natural 
choice if we are to investigate the case of periodic bound- 
ary conditions in SC (a rather natural choice). Thus, any 
relevant physical quantity, /, has the form, 



f(X,Y,Z,T) = f Mm (T) 



JkX+ilY+imZ 



(21) 



Though disturbances are written here as 3D (for further 
studies) we shall focus only on 2D behavior in this work. 
Consequently, we may consider the dynamics within the 
vorticity-streamfunction formulation which greatly simpli- 
fies the analysis. Also, the m subscript on all Fourier com- 
ponents will be suppressed hereafter. 

We now translate relationships I|1U|) - (|12|) into the lan- 
guage of the SC to find, 



u = -d Y ip, 
and 



(9 X 



(d x - 
qTd Y 



qTd Y 



(22) 
(23) 



cU = 0. (24) 

By we immediately see that f w (T) = £ H (T = 0) = 
£ ke is a time invariant quantity. It is, therefore, straight- 
forward to demonstrate that introduction of the solution 
in the form (J2J into (|22I23I) gives, 

1 



u m( T ) 



qTt) 1 
il 



Cm 



(k + qT£) 2 + t 
i{k + qTt) 
(k + qT£) 2 + £ 2 ' 



(25) 
(26) 
(27) 
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The above solutions are the same as those derived by 
Chagelishvili et al. I|2003|l and they contain the possibility 
of TG for some initial conditions. Specifically, a maximum 
in the amplitude of u kl exists when T — T max = —k/(q£). 
This maximum is achieved for T > only for combinations 
of k and t in which k£ < and we will refer to these modes 
as transiently growing or leading. On the other hand, those 
modes in which k£ > always show decay for T > and 
we refer to these as trailing. 

Similarly, a maximum in v ke occurs at T = T max + l/q. 
Nevertheless all linear solutions asymptotically die away 
as T — > oo, where in particular, 



Individual Fourier Mode Linear Energy Evolution 



1 

3.2. Energy 



as, T — > oo. 



(28) 



We define the domain integrated disturbance kinetic en- 
ergy 1 per Fourier mode to be 



E, , — — 
m 2 



+ c.c. 



+ ke i(fcX+ff) + ".) \dXdY, (29) 
which, after using the derived solutions above, reduces to 
\l 12 



E, 



{k + q Tiy + p LxLv -' 



where L x , L y are the periodic scales in the respective X 
and Y directions. 

The transient growth of energy is demonstrated very 
clearly in the expression for E kt : like for the u kt distur- 
bances - a mode will achieve a maximal growth of its en- 
ergy at T = T max and this happens for T > only if 
k/t < 0. 

Let "d kl measure the amplification of the energy of a 
Fourier mode at times T = and T = T max , 



a __ E ke (T max ) _ k^ 

u E k M P + ■ 



(31) 



This definition most clearly shows that the greatest ampli- 
fication occurs for single modes whose ratio k/i is greatest. 
Typical behavior of this transient growth is displayed, for 
several values of k and £, in Fig. ^ 

Letting Ed designate the domain integrated energy of 
the disturbances, it is a simple matter to write it as the 
sum of the individual mode disturbance energies E ke , 



E ri 



v 2 ) dXdY = ^E 



(32) 



Note that this general expression of the domain integrated 
disturbance energy is valid in the nonlinear case - the only 
difference is that the (constant) mode amplitude, £ M , ap- 
pearing in l|30|) will be replaced by the nonlinearly time 

1 Since the kinetic energy is the only form of energy in our 
problem, the "kinetic" qualifier is suppressed hereafter. 




10 12 14 16 18 
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Fig. 1. The behavior of the linear mode energy, E ke , with 
respect to time (1301) in the SC frame is plotted for several k 
and £ Fourier modes. Each mode begins with the same energy 
at T — 0. Transient growth is strongest when the value of \k/£\ 
is largest (see the definition of fl^e in Eq. I3H . The times of 
maximal growth (T = T max ) in the TG modes are also labelled. 

varying amplitude £ ht (T) instead (see below) . The domain 
integrated disturbance energy will follow the same pre- 
scription as in H32|) . 

4. 2D Nonlinear Behavior 



(30) 4.1. Equations Solved and Numerical Method 



We consider fully nonlinear two dimensional disturbances 
(tTO with d z = 0, w = and p' = 0) of the fluid with 
periodic boundary conditions in the SC frame. The equa- 
tions of motion for the disturbances are evolved using the 
streamfunction-vorticity formulation discussed in the pre- 
vious sections. The full nonlinear equations in the NSC 
frame are, 



<9 t £ - qxd y £, + d y i>d x ^ 



Re" 1 (d% + fig£) 



(33) 
(34) 



When the transformations defined by are applied, the 
equations in the SC frame are, 

6U + ^(cU + qTd Y ()-d Y t,(d x Tj) + c(Td Y Tl>) = 

Re" 1 [%£ + (d x +qTd Y ) 2 t] , (35) 

dli>+{d x + qTd Y f^ = i, (36) 

where we have added here, for the first time, an explicit 
viscous term to the equations. Specifically, the equations 
contain now a parameter (the Reynolds number) . The pur- 
pose of this inclusion is to investigate how the flow behav- 
ior depends on Re and what influence Re has on the tran- 
sition to the RTG state. For the most part, however, we 
report on inviscid results here, that is, with formally Re 
= oo in the above equations. We note that the numerical 
procedure implemented (sec below) has a certain amount 
of artificial dissipation associated with it so that though 
the Re number may be infinite, it is so only nominally. 
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Doubly periodic boundary conditions are imposed 
upon all the quantities in the SC frame. This is identical 
with the conditions used by Hawley, Gammie & Balbus 
iP§5)l . 

Equations (|35I36|) are solved using standard Fourier 
spectral methods (Canuto, et al. 1988) except for a num- 
ber of details discussed below. Runs were performed 
usually on specific periodic domains (L x ,L y ) = (ir,2ir). 
Because the nonlinearities are quadratic, a 3/2 dealiasing 
rule is imposed (i) at each time step, after the nonlinear 
terms are computed and (ii) prior to each remapping event 
(see below) . Aliasing problems can be especially hazardous 
here because one can inadvertedly transfer power from de- 
caying (trailing) modes into transiently growing (leading) 
modes and create a situation in which there is spuriously 
generated RTG behavior. The dealiasing procedure out- 
lined above satisfactorilly avoids this problem. 

The Fourier resolution typically utilized 256 modes in 
Y and 512 modes in X. Convergence of the numerical 
scheme was verified by comparing the results of represen- 
tative runs with those of an even higher resolution simula- 
tions (e.g. at 512 x 1024). These numerical tests convince 
us that the 256x512 spatial resolution is satisfactory for 
our purposes at this stage. 

The governing PDE's were temporally evolved using a 
modified Crank-Nicholson method described in Appendix 

m 

The domain in the x — X direction nominally lies 
between X = and X = L x . In the NSC frame, the 
background flow is not moving at x = while it moves 
with speed qL x at x = L x . Consequently this means that 
the shear carries all the y points at x — L x back to their 
initial starting position after a time T rm defined to be 



T. 



qL x 



(37) 



However, as noted in Hawley, Gammie & Balbus (1995), 
points, as carried by the shear, are not always periodic in 
x: they are so only after every integer multiple of T rm . 

Thus, following the general prescription provided for 
by Cabot (1996), the computational domain is remapped 
after every T rm period of time according to the following 
prescription: 

1. At T = t = initial vorticity (or streamfunction) data 
are assumed to be known and given in the NSC frame. 
At this time the NSC (x,y,t) and SC (X,Y,T) are 
coincident. 

2. Given this coincidence of both coordinate systems, the 
initial NSC data are used as the seed vorticity (or 
streamfunction) that evolve by (|35|> and l|36|) , the equa- 
tions in the SC, up to time T — T rm . 

3. Because only at T = T rm (and every integer multiple 
thereof ) are solutions exactly periodic in the x direc- 
tion of the NSC frame, at T = T rm the solution in 
the SC frame is mapped back into the NSC frame by 
applying the inverse transformation of l|16|) onto the 
modes in Fourier space according to the procedure in 



Appendix C, section 2. This act is what we will refer 
to as remapping. 
4. We can treat the mapped, exactly periodic, solution in 
the NSC frame as the new initial condition to evolve 
I|35l36ll which, of course, is in the SC frame. Time in 
SC is reset to zero, i.e. T = 0. The procedure repeats 
with item 2. Note that though T is reset, time in the 
NSC still continues onward. For instance, at the n th 
remapping, T = while t = nT rm . 

This procedure simply amounts to taking a solution 
that has been evolved in the SC and re-expressing it in 
terms of the coordinates of a normal (that is, in NSC) ob- 
server. Then the data as seen by the normal observer is 
used again as initial conditions for the next spate of evolu- 
tion in the SC frame. It should be noted that if there were 
infinite resolution at our disposal, the remapping proce- 
dure mathematically effects nothing since it is a simple 
coordinate transformation. The rationale for applying the 
remapping procedure is technical: it involves issues con- 
cerning finite resolution, periodic boundary conditions, 
and the ability of resolving coherent structures for an ex- 
tended period of time in the SC frame. For more details 
we refer the reader to Cabot {1996) and to the discus- 
sion presented in Appendix C. Prior to each reemapping 
act, the state of the solution is padded using the similar 
3/2 dealiasing procedure (see above). Failing to do this 
can introduce an aliasing error in which a decaying (trail- 
ing) Fourier mode can spuriously turn into a transiently 
growing (leading) Fourier mode (cf. end of section 3.1). 
Nonetheless, we have also performed a comparison run 
to see whether or not the recurrent transient growth we 
observe in the nonlinear solutions (see section 4.3) indeed 
persists when remapping is not done. Aside from slight dif- 
ferences in the evolution of the total disturbance energy 
and the differences in the enstrophy (when remapping is 
applied enstrophy decay is expected and observed while 
when no remapping is performed enstrophy conservation 
is predicted and observed) , the observed recurrent dynam- 
ical phenomena persists in the non-remapped case just as 
it does in the simulations where remapping is done. 

We define the parameter e to be a rough measure of 
what we shall call the turbulence intensity of the system. 
In particular we identify it to be the ratio of the domain 
integrated energy contained in fluid disturbances (i.e. the 
energy Ed as defined by Eq. 13 2|) to the domain integrated 
energy contained in the steady shear flow itself. The latter 
is the domain integral defined by 



J U 2 dxdy = ~f q 2 x 2 dxdy = ^q 2 L 3 x L y . (38) 



E shear = ~ I Vdxdy = ~ 



The turbulent intensity is then formally written as, 
E d 6E d 



E s 



q 2 LlL y 



(39) 



All simulations in this work, unless otherwise noted, 
begin with white noise in the vorticity field with a value 
of e always less than 0.01. That is to say, the energy in 
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the disturbances are always less than one percent of the 
total energy in the steady shear. 

4.2. Disturbance energy and enstrophy evolution 

To gain a certain amount of intuition for the underlying 
processes and to monitor the performance of our numer- 
ical technique we consider the evolution of the domain 
integrated energy of the disturbances. It is a straightfor- 
ward matter to show that in the SC the following holds, 

^± = q J uvdXdY - Re~ l J \Vu\ 2 dXdY , (40) 

where the operator V is given as {dx + qTo\}X. + <9y Y 
and the integrals are taken over the periodic domain. 
The result in (|40l) is known in more general terms as 
the Reynolds-Orr relationship. When Re is finite it clearly 
demonstrates the decaying behavior due to viscosity since 
the second term on the RHS of (|40|l is a negative definite 
quantity. In the limit where Re is so large that the viscous 
decay time scale is very long, the evolution of is gov- 
erned by the first term on the RHS of l|4U|l . It says that for 
q > (< 0) there will be a rise in the integrated kinetic en- 
ergy of the domain if there exists a positive (negative) cor- 
relation between the shearwise and streamwise velocities. 
From H40JI it can be shown that the instantaneous relative 
growth rate of the disturbance energy, that is, din E^/dt is 
independent of the amplitude (see Hennigson 1996). Thus 
the growth rate of a finite amplitude disturbance is essen- 
tially given by mechanisms present in the the linearized 
equations. Without linear growth mechanisms there can 
be no growth even in the nonlinear regime. Nonlinear pro- 
cesses can only shift power among the modes and this, as 
we shall explain below, is a key feature in understanding 
our results. 

Additionally we also consider the enstrophy, 



Z = J \i 2 dXdY. 



(41) 



It too is a straightforward matter to show that this quan- 
tity follows 

% = -Re" 1 J 'iVtfdXdY. (42) 

Thus Z always shows decay for finite Re, irrespective of 
the linear shear, and, it is only when Re = oo is Z con- 
served. In the Re = oo results we get, the enstrophy is 
indeed conserved between remapping events but, as is dis- 
cussed in Appendix C, the act of remapping always re- 
moves enstrophy. The majority of this enstrophy is re- 
moved from very large wavenumber shearwise directed 
Fourier modes and since these contain very little of the 
energy - the overall effect on the global energetics are min- 
imal. 

4.3. Results: RTG and Coherent Vortices 

We have conducted a number of runs for various values of 
the Reynolds number and initial conditions (the spectral 



structure of the disturbance was always the same, but we 
have varied the initial intensity, s(0)). Generally speak- 
ing, we have found RTG behaviour for Reynolds numbers 
above some threshold. In some cases (not high enough 
Reynolds number or too small disturbance strength) the 
RTG was transient, ultimately decaying after a finite time 
(whose size depended on the two parameters mentioned 
above). Inviscid runs (Re=oo nominally) seeded with ini- 
tial e = 0.01 displayed RTG displaying only an insignif- 
icant decay in times of the order of the run duration. 
Consider first Fig. [3 which shows the results of a Re = 
50,000 run, demonstrating the persistence of RTG phe- 
nomenon in the time segment between ~ 70 and ~ 170. 
Long after the linear evolution of the initial conditions 
have died away, the nonlinear behavior shows quasi-steady 
long-lived activity, characterstic of what we have called 
RTG. Guided by the consequences of the Reynolds-Orr 
relationship (see the discussion following Eq. 140ft we have 
performed numerical experiments which demonstrate that 
the peaks in the time evolution are actually driven by lin- 
ear TG. We are led to this proposition by the results of the 
following experiments. At several time points during the 
late evolution of the flow, which are found (in the fully 
nonlinear simulation) to lie just before a growing peak, 
the state of the full flow is read and used as initial con- 
ditions for a parallel numerical evaluation of the linear 
behavior. This exercise reveals that the linear and nonlin- 
ear disturbance energy fluctuations of the full system, at 
least locally in time (between 5 and 10 time units), are 
nearly the same. Eventually the linear flow decays (as it 
should) while the nonlinear flow eventually demonstrates 
TG again. This quality is demonstrated in Fig. El for seven 
separate starting points. 

As stated at the outset, we find that the repeated lin- 
ear TG fluctuations in the total disturbance energy is a 
generic feature of the flows that have been numerically 
investigated by us. The only limitation appears to be 
whether or not there is sufficient power in those set of 
modes of the system that are transiently growing. For in- 
stance, we ran a simulation (not shown) in which all the 
disturbance energy was distributed amongst modes that 
cannot transiently grow. The resulting nonlinear dynamics 
showed fast decay exactly in line with linear theory 

It is also important to note already here that these sys- 
tems, in which sustained RTG is shown, posses an inter- 
esting spatial characteristics. The random initial vorticity 
disturbance develops into collections of coherent long-lived 
vortices which, once developed, translate in the stream- 
wise direction with the speed of the local shear profile. 
It also appears that the number of emerging structures, 
showing temporal persistence, is a function of the box size 
in the shearwise (i.e. disk radial) direction only. This be- 
havior is further detailed below. 

We turn now to a more detailed description of fully 
inviscid (save only for numerical viscosity effect because 
of finite resolution and remapping) runs. As said before, 
these runs show RTG persisting for the full duration of 
our calculations. Three parallel evolution calculations of 
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Fig. 2. Full nonlinear evolution of disturbance energy at Re 
= 50000. Also plotted is the linear evolution of disturbances 
whose initial condition is the fluid state of the nonlinear solu- 
tion at the specified times t = 70.0, 88.0, 102.6, 109.3, 112.0, 
133.3, 146.7. These linear start times are indicated on the plots 
by diamonds. Though the full solution show strong nonlinear ly, 
the linear evolution shows that the rises and falls of the total 
energy appear to be dominated by the linear TG process. 

three nearly identical flows at Re = oo were performed. 
We refer to each of these "runs" as "a" , "b" and "c" and 
their parameters, describing only differing periodic box 
scales, are the following: 

Run a : L x — tt, L y = 2tt, 
Run b : L x = 2tt, L y = 2tt, 
Run c : L x = tt, L y = Air. 

Runs b and c are seeded with the initial data of Run a 
with some additional very weak perturbations. 

Fig. El depicts the evolution of the quantity e for Run 
a. As before we see that the system demonstrates RTG 
phenomenon and the energetics of the system remain ac- 
tive far after the linear evolution (not shown) predicts de- 
cay. After the initial spectrum of decaying modes has died 
away (T < 50) , the system's energy reaches a quasi-steady 
long-term sustained level. The same temporal features de- 
scribed here are observed with Runs b and c. 

By comparison, we have also computed the evolution 
of a viscous version of Run a with Re = 50,000. This run 
demonstrated significant decay of the energy already by 
T = 350 (which was the point where the simulation was 



stopped). The mix between decay and RTG activity is 
clearly evident in both runs here and falls in line with 
implications of the Reynolds-Orr relationship (|4()|) . 

We now turn to the description of the spatial features 
of the results. By t = 100 or so, each of the three runs 
settle down into a quasi-steady state described as a collec- 
tion of coherent vortical structures that propagate along 
with the imposed streamwise flow. These structures ap- 
pear (Fig. |2| & to be long-lived in that they neither 
transiently decay nor do they become swallowed up by 
other vortices (see below) for the duration of the simula- 
tions (up to nearly t ~ 800). In this steady state Runs a 
and c support 3-5 coherent vortices while Run b supports 
about 7-8 coherent vortices. 

Fig.0]is representative of the long time spatial behav- 
ior in all the numerical simulations we have performed. 
It depicts the result of prolific vortex-vortex consump- 
tion leading to the production a single intermediate sized 
coherent structure. All simulations show that while the 
quasi-steady state is achieved, small to medium sized vor- 
tices having emerged at closely neighboring values of x 
eventually consume each other because the background 
shear invariably advects the once-separated structures 
into each others vicinity. Many structured vortices may 
survive into the quasi-steady phase of the evolution but 
only one vortex will be associated with any streamline. 
The vortex merging feature observed here is not unex- 
pected since this sort of thing has been observed in other 
2D numerical investigations and has been discussed at 
length by Mc Williams, (|1984[) . Marcus !|1993[> and more 
recently by Lin et al. ( 2003J) . The comparison to these re- 
sults (and thus testing our scheme) was our motivation in 
running the above three runs, differing only in the com- 
putational domains (see also below). 

There are two main difference in the simulations be- 
tween the finite and infinite Re cases. First, as shown 
above, all finite Re simulations show decay of their activ- 
ity even though the interim dynamics are of the RTG type 
seen in the inviscid cases. In these latter cases (Re=oo) the 
decay times in the total energy are essentially infinite since 
no decay was measurable for the duration of the runs. In 
these cases, any decay can be caused only by the numeri- 
cal dissipation attributed to the remapping procedure (see 
discussion in Appendix C). Second, vortices appearing in 
the finite Re situations are less centrally concentrated than 
their counterparts in the Re = oo cases. Fig. shows this 
differing quality between two versions of Run a comparing 
Re = oo and Re = 50, 000. 

In summary, the results show that (a) in the compar- 
ison runs between domains L y — An and L y — 2n (with 
L x = n in both) , roughly the same number of vortices sur- 
vive and (b) in the comparison runs between the domains 
L x — n and L x — 2tt (with L y = 2tt in both) the num- 
ber of surviving vortices is nearly double in the L x = 2tt 
run. These findings offer qualitative support to the hy- 
pothesis, as suggested by P. Marcus (in conversation with 
one of the authors), that the number of emerging coherent 
structures is only a function of the size of the domain in 
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Fig. 3. The quantity e for Run a: white noise initial condi- 
tions for the vorticity, L x — 2n, L y — 2n. A viscous (Re = 
50,000) and inviscid flow, (Re = oo) are depicted for com- 
parsion. Aside from the initial transient readjustment phase 
(T < 50), the inviscid run demonstrates no decay over the 
course of its duration. 
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Fig. 4. The vorticity distribution at t — 320.0 as viewed by a 
non-sheared observer for Runs a,b,c at Re = 50, 000. In each 
figure a quasi steady state has been reached by the solutions. In 
(a) between 3-5 vortices survive while in (b) 7-8 vortices persist. 
In (c) about 3 vortices survive. The results are suggestive of the 
hyptothesis that the number of surviving vortices is a function 
of the radial (x direction) domain size and that for any given 
streamline there is only one coherent vortex associated with it. 

the shearwise (disk radial) direction and that, further, no 
more than one sizeable long-lived vortex can occupy the 
same radial neighborhood or zone. Of course, the gener- 
ality and veracity of this assertion needs to be evaluated 
with further systematic studies. 

5. Discussion and Summary 

5.1. The shearing box approximation, its limits and 
properties 

We have presented semi-rigorous scaling arguments in or- 
der to derive appropriate equations for describing small- 
scale dynamics in the midplane vicinity of rotationally 
supported thin barotropic disks. The procedure admits 
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Fig. 5. Same as Fig.0]except that Re = oo. The concentration 
of vorticity is tighter than in the viscous case. 



Fig. 6. Surface plot of the vorticity at long times. Comparison 
between Re = 50,000 (top panel) and Re = oo (lower panel). 
The concentration and amplitude of the vorticity in the Re = 
oo case is tighter than in the viscous case. 

two limits, the SSB and the LSB, whose choice depends 
on the length scales of interest. In this work we work with 
the SSB equations, because in this first study we limit our 
interest to dynamical lengths which are small compared to 
the vertical scale-height of the disk. In this way we obtain 
a simpler problem, avoiding the complications resulting 
from the inclusion of acoustic modes. The SSB equations 
are essentially that of incompressible fluid flow in a ho- 
mogenously shearing environment with constant Coriolis 
parameter. In fact, in both the LSB and SSB scaling limits 
the Coriolis parameter is a local constant and the shear is 
homogenous. However, by using a shearing box approxi- 
mation (which utilizes the asymptotic scaling constants e 
and 8 - see Appendix A), a physical scale in the disk is 
set. It appears thus that the argument of Balbus (2003) 
(regarding the above mentioned Longaretti's speculation 
on the numerical resolution needed to see instability in 
rotating flows) is not well founded. 

The asymptotic equations describing the two shear- 
ing box scalings, as derived in the Appendix A, relate 
to those used elsewhere in the literature. The SSB equa- 
tions, in their linearized form, are the same as those used 
by Chagelishvili et al. (2003 ). Identical equations, i.e. de- 
scribing the SSB with no a priori density fluctuations, are 
also equivalent to the ones usually used in the fluid dy- 
namics community to investigate inviscid incompressible 
rotating plane Couette flow (Nagata, [1986). In their full 
nonlinear form, the equations for the SSB are equivalent 
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to the equations cited by Longaretti (20(H| where the den- 
sity disturbances are everywhere zero. 

The equations describing the LSB admit acoustic dis- 
turbances unlike the SSB limit. In their linearized form, 
the equations in the LSB are only similar to the equa- 
tions used by Tevzadze et al. ((^003) where, in their work, 
steady state pressure gradients are assumed to be uni- 
form rather than having the positional dependence they 
are understood to have for barotropic pressure profiles, cf. 
I|A.20|) and (|A.21I) . The fully nonlinear equations of the 
LSB are similar to those used by BHSW except that in 
their treatment the steady density and pressure profiles 
are considered to be uniform (avoiding vertical structure 
effects). This seems to be, thus, a kind of an intermediate 
limit between our SSB and LSB. 

5.2. TG: its linear decay and its nonlinear reccurrence 

Both the finite and infinite Re number simulations de- 
scribed in this work show that long after linear theory pre- 
dicts the decay of all disturbances, the total disturbance 
energy exhibits repeated episodes of transient growth. 
This is established by considering the parallel linear evo- 
lution of the nonlinear state at various points in time. The 
condition of the full nonlinear flow is used as the initial 
condition for the linear evolution and we find that the lin- 
ear and nonlinear development closely shadow each other 
for at least 5 time units. Where there is linear TG the 
full flow also demonstrates TG. Fig. qualitatively indi- 
cates that the temporal pattern in the RTG is complex 
and appears to be chaotic. 

Whether a similar result can be expected in 3D, de- 
spite BHSW repeated claims that this is impossible, must 
be decided by 3D simulations. Recent simulations of 3D 
plane Couette flows (without rotation) at relatively low 
Re-numbers by Schmiegel and Eckhardt l|1997(l show sim- 
ilar (to the one described here) behavior of the integrated 
domain energies. While eventual decay is observed in most 
cases, the fluid response is chaotic in the interim. They 
further show that simulations with the same total initial 
energy and Re numbers but with differing initial condi- 
tions can exhibit vastly differing decay timescales - and in 
some instances there is no dissipation over the duration 
of a simulation (Eckhardt, Marzinzik & Schmiegel 1998). 
Though not reported here in detail, we also observe sim- 
ilar sensitivity to initial conditions for simulations with 
finite Re numbers in the way described in their work. A 
detailed exposition of these results will be reported in our 
future work. 

Though the detailed nonlinear feedback mechanism 
has not yet been identified it seems apparent to us that 
the linear TG phenomenon is a central player of this dy- 
namical system and is a key element in the dynamics of 
rotationally supported flows like thin Keplerian disks. The 
findings of this work also lend support to the conjectures 
put forth by Baggett, et al. l(T555|l . Waleffe ifTMTjl . and 
more recently Grossmann f2QQHI): that subcritical flows, 



in the sense that there are no linearly unstable processes, 
are governed by an interplay between linear TG and some 
sort of nonlinear feedback of energy back onto TG modes. 

5.3. Coherent structures 

Our results have some bearing on certain fundamen- 
tal questions pertaining to 2D turbulence. In the work 
of McWilliams (1984), in which simulations of 2D non- 
globally sheared geostrophic flows were performed, it was 
shown that vorticity concentrations at intermediate scales 
exhibit long lifetimes: i.e. times that are long compared 
to their eddy-turnaround times and with sizes larger than 
the dissipation scales. This result inspired the hypothesis 
that it is a typical feature for two-dimensional geostrophic 
flows to sustain these types coherent vortices and, further, 
for them to retain their independence and identity unless 
there is a close (and usually destructive) encounter with 
another vortex. Additional support for this view was re- 
ported by Bracco et al. (|2000|1 using very high resolution 
2D simulations of the same problem. The appearance and 
persistence of vortices at intermediate scales of the com- 
putational domain in our work qualitatively falls in line 
with the suggestion that the persistence of intermediate 
scale vortices are a generic feature of these sorts of flows. 

Furthermore, the appearance of widespread merging 
of smaller vortices leading to the formation of interme- 
diate scale quasi-steady structures is also consistent with 
what has been observed in the previously mentioned non- 
globally sheared simulations by McWilliams and Bracco 
et al.; and also those of the 3D shearing simulations of 
Lin et al. (2003). What the results of this work add to 
the emerging picture of such flows is that with a glob- 
ally imposed shear, intermediate sized vortices not only 
may appear and persist for a long time but they also tend 
to be solitary. In other words, multiple vortices may ini- 
tially emerge but once the quasi-steady state is reached 
and most of the vortex- vortex merging/destruction events 
have taken place, a given point in the shearwise direction 
x (irrespective of any streamwise position y) will have only 
one sizable coherent vortex associated with it . This pat- 
tern is consistent with the speculations made to this end 
by P. Marcus (private communication). 

We find it somewhat encouraging that the number of 
emerged quasi-steady vortices increases when the domain 
size in the radial direction is increased. It suggests strongly 
that the emergence of coherent structures in such flows is 
an intrinsic phenomenon and it is not, say for instance, a 
by-product of the boundary conditions used. By this we 
mean to say that if at the end of our simulations we were 
to find that only one or two large scale vortices survived 
into the quasi steady state stage irrespective of the size 
of the computational domain in the shearwise direction, 
then we would probably feel that the pattern of developing 
and sustained coherence would not be indicative of phe- 
nomena intrinsic to these flows but rather an artifact of 
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the boundary conditions employed (those here being SB 
periodicity). 

5.4. Summary 

Though the numerical results reported in this work are 
exclusively of two dimensional nonlinear simulations, we 
think that they may contribute to the understanding of 
subcritical transition in more general shear flows as well 
(including rotating ones). We have found here support to 
the idea that excepting the special circumstance, where 
initial conditions are set up via conspiracy (in the sense 
of the phrase as used by Goldreich & Lvnden-Bell. 11965)) 
to have no spectral power in TG modes, these 2D flows 
are likely to be dynamically active and to exhibit and 
maintain robust vortices. Indeed, in viscous 2D flows such 
activity must invariably be limited to finite times, after 
which such activity and the associated vortices will decay. 
Nonetheless, accretion disk Reynolds numbers are enor- 
mous - much larger than those arising from numerical 
roundoff in our calculations. As such, we may expect cor- 
respondingly enormous decay times (already larger than 
the several hundred box orbit times of our simulations) 
during which the occurrence of new perturbations can not 
be excluded. 

It is thus reasonable to still seriously consider the hy- 
pothesis that thin Keplerian disks may be hydrodynami- 
cally active, where the activity we envisage is of the RTG 
kind. It still remains to be seen, of course, whether or not 
three dimensionality and its introduction of stratification 
and Coriolis effects in rotating flows significantly alters 
the implications of our results. 

Based on the results of their simulations, BHSW 
strongly advocate just that, but our work indicates that 
this issue deserves additional consideration. In particu- 
lar, we are now engaged in extending our simulations to 
3D, where different numerical schemes (spectral methods) 
in high spatial resolution and long integration times are 
employed. The effects of stratification and compressibility 
are also being addressed. Since we have found here that 
the initial conditions may be crucial in determining the 
persistence time of RTG, we intend to experiment with a 
variety of these and ascertain optimal perturbations. The 
results of these simulations may decide the important and 
still controversial issue if complex dynamical turbulence- 
like activity may be driven by purely hydrodynamical pro- 
cesses of the RTG kind in rotationally supported thin cir- 
cumstellar disks. 
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Appendix A: Derivation of the Shearing Box (SB) approximation equations 

The shearing box (or sheet) (SB) approximation has been used in a variety of astrophysical applications. It appears 
that the studies of Hill (1878) on lunar dynamics have provided the basic ideas for its development. The approximation 
has more recently been utilized by Goldreich & Lynden-Bell (1965) and Toomre (1981) in the context of galactic disks, 
Wisdom & Tremaine (1988) in the study of planetary rings, Goodman & Ryu (1992)for accretion disks and by Hawley, 
Gammie & Balbus (1995) in the study of magnetized accretion disks. It is within this approximation that most of the 
works mentioned in the Introduction have been done and it is used in this and our subsequent works. It is important, 
in our view, that all the approximations made in deriving the SB equations in their various limits be precisely spelled 
out and clarified and this is the purpose of this Appendix. 



A.l. Global equations in a rotating frame, cylindrical coordinates 

We start with the hydrodynamical equations, describing an ideal fluid (actually an inviscid and isentropic flow), in 
cylindrical coordinates (unit vectors r, 0, z) and in a frame of reference rotating with a given angular velocity Qqz 
(see below). The only body force acting on the fluid is the gravitational attraction of a point mass M, located at the 
origin. The equation of state is assumed to be one of a perfect gas with constant adiabatic coefficient 7. The equations 
of mass momentum and energy conservation are then 

d tP + V- (pw) = 0, (A.l) 

d t w + (w • V)w + 2tt (w r ip - w v r) = — -VP - rQ2 r - V6, (A.2) 

p 

d t P + (w • V)P + 7PV • w = 0, (A.3) 

where p, w and P are the density, velocity and pressure fields respectively and 0(r, z) is the gravitational potential 
which is explicitly given by 0(r, z) = GM(r 2 + z 2 )" 1 ^ 2 . 

Before making the crucial step in the derivation of the SB equations (i.e. considering a small region around a point 
and deriving the local equations valid approximately in this region) we note that for any given rotation law Q = f2(r), 
there exists a steady axisymmetric solution, with the pressure satisfying a global barotropic relation (i.e. being a well 
defined function of p - P — P(p)) and the velocity being composed of the rotation only, that is w = r [Q(r) — flo] tp, 
as expressed in the rotating frame (see e.g. Tassoul & Tassoul 1978). Denoting by Pb{r, z) the density in this solution 
for a given rotation law f2(r) we can express the gravitational force per unit mass as follows 

= -— VP fc + rn 2 r, (A.4) 
Pb 

where P,(r, z) is obtained from p^, by the barotropic relation. Relation l|A.4|) allows one to substitute for the gravity 
term in (| A. 2|> and we shall do it below. It is important to notice that if the rotation law fl(r) is close to the Keplerian 
one and the disk the fluid is cold enough (that is the vertical scale height is small), the steady solution is an essentially 
rotationally supported disk. 



A.2. Local nondimensional box equations in a rotating frame, Cartesian coordinates 

We concentrate now on a small region around a point rg - (r = r*o, 92 = <po, z = 0) - , that is, some typical point in 
the midplane of our disk. Let the region be of size A, so that it (henceforth referred to as the box) is defined by 

A A A AAA 

r --<r<r + -; <p - — < v < ^ + — ; - - < z < -, 

where 5 = — < 1. 

ro 

The smallness of S allows the employment of Cartesian coordinates, since locally the curvature in the azimuthal 
direction is only slight and, as we shall see below, drops out at lowest order in S. Thus we can transform the coordinates 

x = r-r a , y = r (ip - ip ); (A. 5) 



and effect this transformation on equations HA.1IA.3|) , with the gravitational term in (|A.2(I substituted from (|A.4I . 

Before writing out these equations in detail we nondimcnsionalize them, because only then we will be able to 
systematically neglect terms which are of higher order in the small parameter <5. Thus x, y, z are scaled by A = 5ro, t 
by the Kepler time at r , t k = f^^o) = y/r%/(GM), £1 by H^ro), p by p Q = p b (r , tp ,0), P by the corresponding 
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d t w z + (w-V)w z = -(I) (-d z P--d z P b J, (A.9) 



Po and the velocities by vo = Six A = <5f2x(ro)ro. Note that the velocity scale is related to the typical local sound 
speed scale, c s q = \/Po/po by the relation 

S 

vo = -c s0 , 
e 

where e = c s o/[^_R"(^o) r o] an d thus H = era is the local vertical scale- height of the disk. In cold disks, by definition 
e<l and thus they are thin. We assume throughout this discussion that disks are cold. 

We do not specify for the moment the order of magnitude of the ratio of the small parameters 8 and e. Their 
ordering will correspond, as we shall see later, to different physical regimes of the local problem. Thus keeping the two 
small parameters as they are, we obtain the following set of nondimensional equations, which are valid in the small 
box around the point ro: 

d t p + V • (pw) - 0, (A.6) 

d t w x + (w • VK - 2O w H - 2gfi§a; = - (^) 2 {^^P - y b dxPb ) ' ( AJ ) 

/e\ 2 1 

d t w y + (w • V)w y + 2Q w x = - (jj -d v P, (A.8) 

2 fl 1 
' -d t P- — 
,P Pb 

d t P + (w • V)P + 7PV • w = 0, (A.10) 

where we have kept only terms in lowest order in 5 and therefore neglected the curvature terms, so that the above 
equations (that is, the vectors and derivatives) are actually expressed in local Cartesian coordinates, defined in (|A.5|) . 
The net radial body force in the rotating frame r(fl 2 — fl^) has been expanded around ro up to order S, where fio 
the rotational velocity of the frame is actually the value of the rotational law angular velocity at ro, that is, fi(ro) 
expressed in units of the Keplerian rotational velocity at that point. For Keplerian rotation law we obviously have 
fio = 1. The parameter q is related to the Oort constant A: 

In the terms on the right hand side of the momentum equations we have the differences between the pressure force 
components and the corresponding components in the basic steady barotropic state, referred to above. That barotropic 
state is clearly axially symmetric and thus d y P\, — d y pb — 0, but in the radial and vertical directions we can not expect 
a priori the uniformity the density and pressure. 

Note, however, that because the ratio (e/<5) 2 multiplies the pressure force terms in question, their size will have to 
be such as to allow for proper balancing of the equations (see below). 

A. 3. Removal of steady state linear shear - the shearing box (SB) equations for the perturbations 

We notice now that w = 2Axy, P — Pb and p — pb is a steady solution of equations (|A.6IA.10J| and therefore we 
consider this solution as our basic flow and perturb it by adding functions of space and time (denoted by primed 
quantities). Thus we write 

w = 2Axy + u', (A.12) 
P = Pb + P ', (A.13) 
P = P b + P', (A.14) 

Substitution of the above form of the functions in the equations gives rise to the removal of the linear shear steady 
solution and the following equations for the perturbations u', p' and P' are obtained: 

dtp' + 2Axd v p' + Pb V ■ u' + u' z d zPb + u' x d xPb + V • (pV) = 0, (A.15) 

cK + 2A xdy< - 2o 0U ; + k . v K - - (|) a - ^ ) , (A.16) 

d t u' y + 2Axd y u' y + 2(A + Q )u' x + (u' ■ V)< = - (ff ) ' ( A ' 17 ) 

d t u' z + 2Axd v u' z + (u' • V)< = - Q) 2 - ) ■ (A.18) 

\oJ \Pb + p' Pl + PbP J 

d t P' + 2Axd y P' + u' x d x P b + u' z d z P b + 7AV • u' + u' ■ VP' + 7 P'V • u' = 0. (A.19) 
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where wc have retained the full nonlinear terms, that is, did not expand in the perturbations. 

This set of equations for the perturbations is quite general and it includes all the different cases which were 
used in the literature. The advantage of this formulation is in its consistency and physical transparency as to which 
approximations are made (see below). 



A. 3.1. Small shearing box (SSB) equations (5 <^ e <C 1) 

When we chose 8 -C e, the shearing box size is much smaller that the disk's thickness (the unperturbed pressure 
vertical scale height). This choice implies that the pressure force terms in equations ljA.16IA.18p are multiplied by a 
very large factor and therefore if we want that in the prevailing dynamics these terms be of the order of the inertial 
terms, they themselves must be of the order of (8/e) 2 . 

Consider first the terms containing the gradient of the steady barotropic pressure P,. The z derivative of Pb (and 
thus of pb as well) is in our units clearly of the order (8/e) 2 , because by nondimensionalizing equation (|A.4|I one gets 
that in the box 



~j {1 + Sx){n 2 (x)~ [(l + 8x) 2 +8 2 z 2 ] 3/2 }, (A.20) 
l -d z P b =-(-) 2 z [(1 + Sx) 2 + 8 2 z 2 } " 3/2 , (A.21) 



Pb \e- 



Pb \e 

where f2(x) is the rotational nondimensional rotational velocity in the box. 

The order of magnitude of j^d x Pb depends on the rotation law. If the rotation law is exactly Keplerian, we have 
fl 2 (x) = (1 + Sx)~ 3 and —d x P b = O (8 3 /e 2 ) and —d z P b — -z8 2 /e 2 + O (S 3 /e 2 ). More generally speaking, in order to 

Pb ^ ' Pb * ' 

have j^d x Pb = O ((5 3 /e 2 ), the basic rotation state does not have to be exactly Keplerian, but just close enough to it. 
In other words, order 1 departures from a Keplerian rotation profile are needed to see significant radial gradients of 
the steady state in this approximate limit. 

Thus, assuming that the rotation law is close enough to the Keplerian (some special choices appear in the body of 
the paper) , and for the sake of economy of the notation, we make the replacements 

d l P b (-) d t P b (A.22) 



and 



2 



d lPb -»• I - I d lPb (A.23) 



for i — x, z, remembering that the units of the barotropic profile spatial derivatives have been scaled accordingly. 

Examine next the terms containing the perturbations P' in equations HA.16IA.18|) . If these are not sufficiently 
small (that is of order 8 2 / e 2 ) we may formally write 

P ' = P o+ (£) P+- ( A - 24 ) 

With this expansion (and using the above scaling of the barotropic terms) , equations ljA.16IA.18p give in lowest order 

VP^ = (A.25) 

that the O (1) part of the pressure perturbation is uniform in the box. Using this in equation (jA.19jl gives in lowest 
order 

d t P^ = -jV-[(P b + P^u'}. (A.26) 

Integration of this equation over the volume of the box shows that the time variation of Pq depends on the boundary 
conditions. For both periodic or zero velocity boundary conditions (which are usually employed) we get that Pq is 
also constant in time. Thus this constant pressure shift may be formally absorbed into P, and consequently changing 
nothing. Thus we may now formally write 

P' = (J\ p, (A.27) 

in which p is now understood to be an order one quantity. This actually amounts to redefining the units of the pressure 
perturbation. 



17 



Substitution of these scalings into (|A.15IA.19|I and dropping terms of O ((5/e) 2 ) and smaller gives that the energy 
equation (|A.19(1 simplifies to 

V • u' = 0, (A.28a) 

that is, we are dealing in this approximation with an incompressible flow. This is physically consistent with the fact that 
the relative pressure perturbation was assumed to be very small. We have an incompressible flow with the density of 
any fluid element conserved by the flow and the the small pressure perturbations describe the dynamics of non-acoustic 
(that is solenoidal, or vortical) modes only. 
The three momentum equations become: 

d t u' x + 2Axd y u' x - 2Q u' + (u' • V)< = - ( — T~, dxP - -5-? 7 d x P b ] (A.28b) 

\Pb + P' Pt + PbP ) 

d t u' y + 2Axd y u' y + 2(n + A)u' x + (u' • V)u' y = - ( - ; d y p) , (A.28c) 

d t u' z + 2Axd y u' z + (u' • V)u' z = - (—^— d zP - d z P b ) (A.28d) 

\Pb + p' Pt+PbP J 

Finally, the equation for the density perturbation p' results from (|A. 15|) and, in this order, reads 

dtp' + 2Axd y p' + u' • V// = 0. (A.28e) 

The set (IA.28alA.28ell . supplemented by the appropriate boundary conditions (usually chosen to be periodic), describes 
the full nonlinear problem in the SSB approximation, that is, for S <C e. These equations with p' = throughout is 
equivalent to the system studied by Longaretti 2002 and is often times referred to and studied as rotating plane 
Couette flow (Nagata, 1986|. The linear problem contains the one studied by Chagelishvili et al. (2002). 

A. 3. 2. Large shearing box (LSB) equations (5 ~ e <C 1) 

We chose now 6 = e the size of the shearing box is of the order of the vertical thickness of the disk. In this case the 
relative size of pressure perturbations is of the order of the relative density perturbations (to balance the relevant 
terms) and the vertical and horizontal derivatives of the barotropic pressure are necessarily retained in all of the 
equations as well. When the flow is nearly Keplerian, only the vertical derivatives make it in at the lowest order as it 
was argued in the previous section. We get equations identical to the full set ljA.15IA.19f) but with S = e and where, 
in practice, only the first order terms of the basic state barotropic profiles l|A.20IA.21|) survive to this lowest order. 
Unlike the SSB, in this limit acoustic waves will be included in the equations as the flow is not incompressible. We 
shall treat this system in future works. 

It should be noted that Tevzadze et. al (2003) used a simplified version of the linear limit of these equations. The 
system studied by Balbus, Hawley and collaborators (for instance Eqs. 3.5a-b in Balbus, Hawley and Stone, 1996) is 
a special case of this problem with the background barotropic density and pressure profiles assumed to be uniform. 

Appendix B: Details of the nonlinear 2D numerical method 

We begin with and (23). We assume all quantities to be Fourier-Galerkin expanded in both the X and Y directions. 
For example we have for the vorticity, 

k,e 

we similarly expand in this way for the stream function ■0. Note that each In the Fourier basis the equations are in 
which the Jacobian J is 

J(il>,t)=d Y il>(d x t + <frd Y t)-d{ Y (d x il> + <frd Y 1>). (B.2a) 

Like £ fcf , J k( references the k,£ Fourier component of the Jacobian. The time stepping is done via a modified Crank- 
Nicholson method outlined below. We write the above equation in the more generalized form, 

f* = -'u+£ M (T)£ M) (B.3) 

where the linear operator C M acting on Fourier component k has an explicit T dependence. The evolution scheme 
time centers the evaluation of the nonlinear term J ki and uses a time averaged evaluation of the linear operator. If 
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superscript n refers to the n th time step and if St refers to the time difference between successive time steps then the 
temporal discretization of l|B.3(l becomes, 



C^-C" 1 CUT + St^+C^T-St)^- 1 

25t = 2 J - (T) ' (R4) 

reshuffling terms gives, 

tn+i = 1 + StC ke (T - St) x 2St 

I - StC kl (T + St) ^ l-St£ kt {T + St) kA ) ' { ' 

We also remind the reader that J™ is the Jacobian evaluated at the time step n. We notice that the terms involving the 
linear operator appearing in the above expressions are the first terms of a Taylor series expansion of their exponential 
forms. That is, 



exp { / C M dT' Ul + StC ki (T) + O (St 2 ) , (B.6) 




we replace all appearances of C kl with the above approximate form. Thus (|B.5|) becomes 

C pT+St ~\ C pT+St ~\ 

C +1 = ex P \ j T st ^dT' \ CJ 1 - 25t exp I ^ C kt dT' \ J£ (T). (B.7) 

Replacing the linear operators as we have here is akin to a Pade Approximation (Bender & Orszag, 1999) except in 
reverse: we have "inferred" a more accurate form of the operator by recognizing that the fractional terms in (|B.5|) are 
the approximations of the exponentials which are the exact solutions to problems with such operators. This procedure 
more accurately represents the linear evolution and it removes the well-known high-wavenumber inaccuracies (though 
stable) associated with Crank-Nicholson methods: in particular, as originating from the denominators of (|B. 



Appendix C: Remapping: its rationale and technical details 

The original call to our attention of the ideas behind remapping is credited to P.S. Marcus who, in several private 
comminications with one of the authors (OMU), sketched out the essential rudiments of the problems (and their 
resolution) as they are discussed herein. To our knowledge, Cabot (|1996(l is the only published work in the literature 
in which the steps describing the remapping procedure is both outlined and discussed with some detail. 

C.l. The rationale 

We have two goals in mind when developing (numerical) nonlinear solutions of (|35|l and (|36|) : 

f. to accurately implement periodicity in the SC frame (as in the programme outlined in Hawley et al.,[lp95.) and, 
2. to resolve spatially coherent structures for as long a time as possible. 

The first goal is automatically achieved by going into the SC frame and by assuming a spatially doubly periodic 
Fourier-Galerkin decomposition of the solution (see Appendix B). 

The trouble lies in the inability of resolving spatial coherence for an arbitrarily long time in the SC frame. To 
appreciate why this a problem in the first place, we must first define what we mean when we speak of spatial coherence. 
For a non-sheared normal observer, i.e. someone in the NSC frame, we understand a coherent structure to be a fluid 
structure that has a more-or-less definite size and shape. By long-lived we imply that the coherence in the structure 
remains constituted as such for a long time, at least in the statistical sense. In the language of this article, it is to say 
that in the NSC frame the coherent structure will have a statistically-steady wavenumber bounded spatial spectral 
profile with power peaking at some finite set of Fourier wavemodes. For the purpose of this argument, the mechanism 
behind the statistically steady maintenence of power at these wavemodes is inessential. For the sake of transparency 
in the subsequent discussion, let us suppose that all the power resides exactly with one wavemode with wavenumber 
ko = (ko,£o) in the NSC frame. 

The problem centers around the temporally limited ability to represent in the SC frame an object that appears 
(statistically) steady in the NSC frame. As time progresses the representation of a coherent structure in the SC frame 
steadily degrades fundamentally because the computational resolution in the SC frame is limited. Below we describe 
in more detail the origins and nature of this problem. 
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To express these ideas in more mathematical terms we begin with a coherent structure whose vortical power is 
concentrated, as we said, in a single Fourier wavemode with wavenumber ko,£o expressed in the NSC frame as, 

&,y,t)~L oeo e ikox+Uoy , (ci) 

in which £(x,y,t) is the representation of the vorticity in the NSC frame and where £, koto is the complex amplitude 
of the vorticity at wavemode ko,£o- Because the dynamics are steady we assume that £,k e is time-constant 2 . Let 
n -1 (t ) eg) symbolically represent the act of transforming from the NSC to the SC frame at the reference time t 
(i.e. the mapping defined by the inverse of Eg. I16fl and let II(t ) ® represent the forward transformation. Letting 
Y, T) represent the vorticity in the SC frame, then the forward and backward coordinate mapping procedure on 
the vorticity is formally the following, 

n- 1 ^)®^!,!/,*) — *£(X,Y,T), Il(t )®t(X,Y,T)^t(x,y,t). 

Specifically speaking, applying II _1 (i ) on the vortical mode in IjC.ljl at some starting reference time t in the NSC 
frame gives, 



^(X, Y, T) = ^eiiko-lTeoWioY = ^ik^X+itoY^ kcff _ _ qT ^ (C 2) 



where k e g is an effective X direction wavenumber whose magnitude eventually grows with T. Since t is the reference 
time for the transformation, k e g and fco are initially identical because T — when the result of the transformation is 
viewed at t = t . 

It is easy to see that the coherent structure, as viewed from the SC frame, needs to be represented by an increasingly 
larger set of Fourier modes in the X direction. But because we are computationally limited to finite resolution (even 
in the SC frame) eventually one will run out Fourier modes in the SC frame that can properly represent the coherent 
structure. We assume that in the SC frame the X domain is represented by a finite set of wavenumbers, ki, lying within 
and bounded by a maximum wavenumber ±fcjv, or, in other words, — fcjv < fcj < kpf, with fcjv > 0. The bounding 
wavenumber, kpf, is set by the number of grid points (say N) in the X direction of the domain. 3 The magnitude of 
fc e ff will exceed fcjv after a time T res given by 

rp k : k N 



res 



qio q\£o\ 



For T > T res the numerical simulation in the SC frame requires there to be a Fourier mode |fc| = |fc e ff| > fc/v in order to 
represent the coherent structure. Because this mode is not there the coherent structure can no longer be represented 
within the physical construction of the computation. We refer to this effect as permanent information-loss which is, 
in short, very bad. By contrast, with infinite spatial resolution at our disposal (i.e. fcjv = oo) this information-loss 
phenomenon would not occur because there would always be a k = fc e ff represented in the SC frame for all T > 0. 

In general a coherent structure will be represented by a number of Fourier modes and not the single one assumed for 
illustration here. Nevertheless the effects just described applies to them as well since the representation of their power 
in the SC frame will drift in wavenumber according to the process just described. Thus, in considering a resolution 
of this dilemma, we must always keep in mind that the coherent structure(s) of interest usually have some compact 
spectral representation in Fourier space. In desiring to avoid major information loss about the coherent structure(s), 
we mean to protect and preserve as much as possible the constituents of its power spectrum while it is represented 
and dynamically evolved in the SC frame. 

One final observation is required before the solution may be apparent. Irrespective of resolution matters, one may 
also, after any passage of time At, apply the coordinate transformation n(t ) back upon the solution evolved in the 
SC frame and recover the original steady coherent structure with wavemode (fco, A)) in the NSC frame. In other words 
the transformation, 

n(i ) ® £(X,Y, At) - Si{x,y,t + At), 

recovers the original coherent structure in the NSC frame. For our particular example we would have 

t(x,y,t + At) = £{x,v,i a ) = C koeo e lkoX+Uoy , t = t + At. (C.3) 

One may then restart the evolution in the SC frame by applying n _1 (t ) upon £(x, y,i ) noting that t is the new 
reference time, which is T rm after the original reference time t . This means that the simulation in the SC frame begins 

2 Even this assumption is too restrictive. The amplitude could have some time dependence in the argument that follows. 
For the argument presented, it is enough that the power remains contained in this single wavemode, or at worst, a handful of 
wavenumber bounded wavemodes. 

3 Naively speaking: if the size of the domain in the X-direction is L x then roughly speaking fcjv = N/2L X . 
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anew with T again set to zero and it is valid again for T > 0. Wc call this remapping. In terms of the example we are 
using for illustration, the single wavemode vorticity again starts off with its power in the wavemode k c g = ko as before. 
The benefit of repeated applications of this procedure (of evolving and remapping) is that coherent structure(s) can 
be indefinetly represented in the SC frame so long as the remapping step is done well before the resolution loss time 
is approached in the SC frame. In our example it means that the single wavemode vorticity can be represented for as 
long as one wishes as long as T rm < T res . Without the act of repeated remappings the coherent structure can only be 
represented up to the point T — T Ies and no longer. 4 

The way out of this information-loss dilemma makes use of the observations in the previous paragraph. Starting 
from some time t , one may perform the numerical computation in the SC frame up until a time T rm , which is 
much less than the time T TCS when significant information-loss sets in. At T rm the solution generated in the SC is 
coordinate mapped into the NSC using II(t ) (i.e. Eq. lTo)l . Once back in the NSC we immediately reapply the inverse 
transformation Il~ 1 (t + T lm ) (i.e. the inverse of Eq. ^J|with the reference time set to t + T rm ) upon the solution 
in the NSC frame to begin the simulation again in the SC frame. The evolution of the solution in the SC frame 
proceeds with T reset to zero again. If there are intrinsic steady coherent structures in the flow, this act of repeated 
remappings always corrals their wave power into a wavenumber range that is always resolved in the SC frame; thereby 
preventing the information-loss just outlined. It is very important, then, to apply this remapping procedure well before 
any significant information about the coherent structure is permanently lost, i.e. make sure to ensure that T rm <C T les . 

We summarize the procedure below: 

(i) Vorticity data in the NSC frame at time t , i.e. y,t ), is transformed into the SC frame via the coordinate 
map: 

(ii) The vorticity in the SC frame is evolved from T = to T = T rm . 

(iii) The solution of the vorticity at time T = T rm in the SC frame is expressed in terms of what it would look like in 
the NSC frame: 

H(t.) ® £(X, Y, T rm ) — » £(x, y, t + T rm ). 

(iv) The reference time is updated from t — * t Q + T rm and the steps repeat again with (i) . 

The choice for T rm is technically arbitrary. However, the natural choice is to use those moments in time when the 
solution appears exactly periodic in the (shearwise) x direction in the NSC frame. This is detailed in points 1-4 in 
section 4.1 of the text. 



C.2. Some technical details 

The solution for £ in the SC frame is represented as in (|B.1|I . The wavenumbers in the X and Y direction respectively 
are taken to be 

fc = ^Ps n t e0,±l,±2,-,±n mo:t , e=^^-, m ( e0,l,2,---m mai . (C.4) 

L X Ly 

Because of the complex conjugate symmetry of the solutions in wavenumber space we need only consider half the 
modes (the positive ones including zero) in the Y direction. Because this is a numerical simulation we are limited 
to a maximum number of modes to be n max , m max in the X and Y directions respectively. Between remappings the 
maximum dynamically resolved wavenumbers are ±k max , £ max given by, 

i 27T71 ma;c 'llTUl n i ax 

fcmax — j i ^raax — j 3 
L x Ly 

and where, of course, k max > 0. The vorticity is evolved from t to time t + T rm whereupon we have 

a* +T rm ) =J2^(t +Trm)e ikx+ieY . (C.5) 

k,t 

At this stage the solution is considered to be periodic in the non-shearing frame. Thus we reexpress (|C.5(I in the 
non-shearing frame by using the transformation according to (|16(l to find, 

^ + Trm) - Y.^ t o +T ^y kX+ieV+ ' lqTrmtX 
k,t 

= E f M (*□ + T rm y (k+Ak > +a y , (c.6) 

kj, 

4 Mathematically speaking, the basis for why this is all acceptable lies in the fact that the transformation t — > t + At, where 
At is a constant, leave invariant the fundamental equations I33I34H . 
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in which 

2ir 

Ak = m e —. (C.7) 

Now because at the remap time the two coordinate systems are coincident we may reexpress the solution l|C.6(l in the 
SC frame, 

k,l 

with the reference time t fcf. I16|l updated by a factor T rm : in other words, at this point t — > t a + T rm . A momentary 
inspection of (|C.8() reveals that at the time of remapping, in which the SC and NSC frames are coincident, wholesale 
shift of power occurs from leading, transiently growing modes towards trailing, decaying modes. 5 For example sup- 
posing our computational domain was such that L x = L y = 2n, if we had a wavemode with k — —2 and I — 3 then, 
according to l|C.8|) , after the remap that wavemode would become k = 1 and I = 3 based on the shift required by l|C.7|) . 
As we see, in the updated SC basis, subseqent each remapping event the reference time to (cf. lltj|) accrues an additional 
factor of T rm . Most critically, if the remap shifts power to be outside the trailing edge maximum wavenumber, +k max , 
that power is considered lost to the system. 

The following illustrations, beginning with the example in the previous paragraph, should help clarify all of this. 
Suppose that k max = 10 and we are considering power in horizontal wavemodes with t = 3. As in the previous example, 
power in wavemode £_ 2 3 before remapping becomes £ t 3 after the remap according to (|C.7IC.8Jl . Similarly £ 5 3 — > £ 8 3 . 
Near the trailing edge limit, for example, £ 9 3 would become £ 12 3 , but because the maximum positive (trailing edge) k 
wavenumber is 10 the power in the mode £ 3 gets shifted out of the computation all together and it is, in this sense, 
totally lost. On the other hand suppose we were considering the wavemode £_ 10 3 : after remapping its power shifts into 
wavemode £-7,3. After remapping, however, the £_ 103 wavemode is initiated with zero power. In fact, for the I — 3 
wavemodes, the modes £ 10 3 , £_ 9 3 , £_ 8 3 are all initiated with zero power. Had we been considering instead wavemodes 
with i = 5 then after the remapping the wavemodes £_ 10 5 , £_ 9 5 , £_ s B , 5 , 5 would all be set to zero. 

It should be clearly evident from this that a large amount of enstrophy will be lost to the system after each 
remapping event. Because we can predict exactly which wavemodes prior to remapping will have its power shifted out 
of the computation, we can (and do) predict exactly how much enstrophy (and, for that matter, energy) will be lost after 
each remap. In fact, the numerical method employed in evolving the equations here is such that between remapping 
events the enstrophy is exactly conserved. However, what is more important is that after a remap power is explicitly 
shifted towards more positive values of k. Additionally after a remap, the most extremely transiently growing modes 
(modes with very large and negative k wavenumbers) become initiated with zero power after remapping. Because this 
remapping procedure, as prescribed here, is executed in the Fourier-Galerkin basis and not in physical space, there 
is no potential for aliasing errors here. In other words, the remapping will not artificially create power in transiently 
growing modes. This means that we can be assured that any power that develops in transiently growing modes do 
so naturally (i.e. are truly physical) via the usual three wave interactions understood for such systems like this 2D 
inviscid incompressible flow. 



Recall from linear theory that for wavemodes with i > 0, transient growth is expected for modes with k < and decay for 
k> 0. 
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